Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(broom.mixed) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(ggeffects) #for effects plots in ggplotjk
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(DHARMa)    #for assessing dispersion etc
library(glmmTMB)    #for glmmTMB
library(performance) #for diagnostic plots
library(see)         #for diagnostic plots
library(ordinal)    #for ordinal models
theme_set(theme_classic())

Scenario

Read in the data

hughes <- read_csv('../data/hughes.csv', trim_ws=TRUE) %>%
  janitor::clean_names()#bleaching data for 2016
## Parsed with column specification:
## cols(
##   REEF = col_character(),
##   HABITAT = col_character(),
##   SECTOR = col_character(),
##   SCORE = col_double()
## )
glimpse(hughes)
## Rows: 3,394
## Columns: 4
## $ reef    <chr> "09-357", "09-357", "09-357", "09-357", "09-357", "09-366", "…
## $ habitat <chr> "C", "C", "F", "F", "U", "C", "L", "C", "C", "L", "C", "C", "…
## $ sector  <chr> "North", "North", "North", "North", "North", "North", "North"…
## $ score   <dbl> 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 2…
reef habitat sector score
09-357 C North 4
09-357 F North 4
09-357 U North 3

Data preparation

hughes <- hughes %>%
  mutate(oscore = factor(score, ordered = TRUE),
         habitat = factor(habitat),
         sector = factor(sector, levels = c('North', 'Central', 'South')),
         reef = factor(reef))

For a categorical variable, it will apply a treatment contrast.

hughes.sum <- hughes %>%
  count(sector, habitat, oscore) %>%
  group_by(sector, habitat) %>%
  mutate(prop = prop.table(n),
         oscore = factor(oscore, levels = rev(levels(oscore))))

Exploratory Data Analysis

ggplot(data = hughes.sum, aes(y = prop, x = habitat)) +
  geom_bar(stat = 'Identity', aes(fill = oscore), color = 'black') + #'Identity' means don't do anything on it, multiply by 1
  facet_grid(~sector) +
  theme_bw() +
  theme(panel.spacing.y = unit(10, 'pt'))

Fit models

hughes.clmm <- ordinal::clmm(oscore ~ habitat*sector + (1|reef),
                             data = hughes)

#partial plot
plot_model(hughes.clmm, type = 'eff',
           terms = c('sector', 'habitat'))

#summary
summary(hughes.clmm)
## Cumulative Link Mixed Model fitted with the Laplace approximation
## 
## formula: oscore ~ habitat * sector + (1 | reef)
## data:    hughes
## 
##  link  threshold nobs logLik   AIC     niter       max.grad cond.H 
##  logit flexible  3394 -3399.83 6831.67 1692(14086) 2.40e-03 6.5e+02
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  reef   (Intercept) 7.812    2.795   
## Number of groups:  reef 809 
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## habitatF                -0.3714     0.1930  -1.924 0.054304 .  
## habitatL                -1.4684     0.1909  -7.693 1.44e-14 ***
## habitatU                 1.5678     0.4732   3.313 0.000924 ***
## sectorCentral           -2.2063     0.3630  -6.078 1.21e-09 ***
## sectorSouth             -7.1817     0.3017 -23.807  < 2e-16 ***
## habitatF:sectorCentral  -0.6649     0.3389  -1.962 0.049745 *  
## habitatL:sectorCentral   0.3329     0.3397   0.980 0.327179    
## habitatU:sectorCentral  -1.0853     0.6595  -1.646 0.099831 .  
## habitatF:sectorSouth    -0.3648     0.2494  -1.463 0.143585    
## habitatL:sectorSouth    -0.4088     0.2541  -1.609 0.107574    
## habitatU:sectorSouth    -0.7203     0.5323  -1.353 0.176024    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 0|1  -9.1435     0.2725 -33.553
## 1|2  -5.8773     0.2291 -25.659
## 2|3  -3.1539     0.1971 -15.998
## 3|4  -0.7224     0.1784  -4.049
exp(-9.1435) #the odds for habitat C for having 0 bleaching compared to having a higher bleaching that that is 0.0001069125%
## [1] 0.0001069125
exp(1.5678) #odds of having higher bleaching score than habitat C in North sector
## [1] 4.796085
1/exp(-7.1817)
## [1] 1315.142
emmeans(hughes.clmm, ~oscore|habitat+sector, mode = 'prob') #gives a weighed average so the sum of each column probability is 1
## habitat = C, sector = North:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      1.07e-04 2.91e-05 Inf  4.98e-05  1.64e-04
##  2      2.69e-03 6.12e-04 Inf  1.49e-03  3.89e-03
##  3      3.81e-02 7.18e-03 Inf  2.41e-02  5.22e-02
##  4      2.86e-01 3.28e-02 Inf  2.22e-01  3.50e-01
##  5      6.73e-01 3.93e-02 Inf  5.96e-01  7.50e-01
## 
## habitat = F, sector = North:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      1.55e-04 4.80e-05 Inf  6.10e-05  2.49e-04
##  2      3.89e-03 1.05e-03 Inf  1.83e-03  5.95e-03
##  3      5.42e-02 1.25e-02 Inf  2.97e-02  7.88e-02
##  4      3.55e-01 4.46e-02 Inf  2.67e-01  4.42e-01
##  5      5.87e-01 5.66e-02 Inf  4.76e-01  6.98e-01
## 
## habitat = L, sector = North:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      4.64e-04 1.39e-04 Inf  1.91e-04  7.36e-04
##  2      1.16e-02 3.01e-03 Inf  5.67e-03  1.74e-02
##  3      1.44e-01 2.90e-02 Inf  8.75e-02  2.01e-01
##  4      5.22e-01 2.59e-02 Inf  4.71e-01  5.73e-01
##  5      3.22e-01 5.16e-02 Inf  2.21e-01  4.23e-01
## 
## habitat = U, sector = North:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      2.23e-05 1.20e-05 Inf -1.20e-06  4.58e-05
##  2      5.62e-04 2.89e-04 Inf -4.40e-06  1.13e-03
##  3      8.24e-03 4.07e-03 Inf  2.70e-04  1.62e-02
##  4      8.31e-02 3.66e-02 Inf  1.13e-02  1.55e-01
##  5      9.08e-01 4.09e-02 Inf  8.28e-01  9.88e-01
## 
## habitat = C, sector = Central:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      9.70e-04 3.47e-04 Inf  2.91e-04  1.65e-03
##  2      2.38e-02 7.70e-03 Inf  8.75e-03  3.89e-02
##  3      2.55e-01 5.64e-02 Inf  1.44e-01  3.65e-01
##  4      5.36e-01 2.33e-02 Inf  4.90e-01  5.82e-01
##  5      1.85e-01 4.81e-02 Inf  9.05e-02  2.79e-01
## 
## habitat = F, sector = Central:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      2.73e-03 1.07e-03 Inf  6.35e-04  4.82e-03
##  2      6.42e-02 2.22e-02 Inf  2.08e-02  1.08e-01
##  3      4.55e-01 6.87e-02 Inf  3.21e-01  5.90e-01
##  4      4.03e-01 6.64e-02 Inf  2.73e-01  5.34e-01
##  5      7.45e-02 2.54e-02 Inf  2.47e-02  1.24e-01
## 
## habitat = L, sector = Central:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      3.01e-03 1.24e-03 Inf  5.81e-04  5.45e-03
##  2      7.04e-02 2.55e-02 Inf  2.04e-02  1.20e-01
##  3      4.73e-01 7.06e-02 Inf  3.35e-01  6.12e-01
##  4      3.85e-01 7.20e-02 Inf  2.44e-01  5.26e-01
##  5      6.79e-02 2.49e-02 Inf  1.92e-02  1.17e-01
## 
## habitat = U, sector = Central:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      5.99e-04 3.32e-04 Inf -5.08e-05  1.25e-03
##  2      1.49e-02 7.86e-03 Inf -5.35e-04  3.03e-02
##  3      1.78e-01 7.42e-02 Inf  3.21e-02  3.23e-01
##  4      5.38e-01 2.74e-02 Inf  4.84e-01  5.92e-01
##  5      2.69e-01 1.04e-01 Inf  6.57e-02  4.72e-01
## 
## habitat = C, sector = South:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      1.23e-01 2.17e-02 Inf  8.08e-02  1.66e-01
##  2      6.63e-01 2.02e-02 Inf  6.24e-01  7.03e-01
##  3      1.96e-01 2.99e-02 Inf  1.37e-01  2.54e-01
##  4      1.59e-02 3.41e-03 Inf  9.26e-03  2.26e-02
##  5      1.56e-03 3.80e-04 Inf  8.18e-04  2.31e-03
## 
## habitat = F, sector = South:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      2.27e-01 3.88e-02 Inf  1.51e-01  3.03e-01
##  2      6.58e-01 2.29e-02 Inf  6.13e-01  7.03e-01
##  3      1.07e-01 2.12e-02 Inf  6.50e-02  1.48e-01
##  4      7.71e-03 1.89e-03 Inf  4.00e-03  1.14e-02
##  5      7.49e-04 2.04e-04 Inf  3.50e-04  1.15e-03
## 
## habitat = L, sector = South:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      4.79e-01 5.59e-02 Inf  3.69e-01  5.88e-01
##  2      4.81e-01 4.82e-02 Inf  3.87e-01  5.76e-01
##  3      3.71e-02 8.58e-03 Inf  2.03e-02  5.40e-02
##  4      2.48e-03 6.57e-04 Inf  1.19e-03  3.77e-03
##  5      2.40e-04 6.94e-05 Inf  1.03e-04  3.76e-04
## 
## habitat = U, sector = South:
##  oscore     prob       SE  df asymp.LCL asymp.UCL
##  1      5.68e-02 1.59e-02 Inf  2.56e-02  8.81e-02
##  2      5.55e-01 5.32e-02 Inf  4.51e-01  6.60e-01
##  3      3.48e-01 5.68e-02 Inf  2.36e-01  4.59e-01
##  4      3.63e-02 1.03e-02 Inf  1.60e-02  5.65e-02
##  5      3.64e-03 1.15e-03 Inf  1.39e-03  5.89e-03
## 
## Confidence level used: 0.95
emmeans(hughes.clmm, ~habitat|sector, mode = 'mean.class') #mean bleaching for each sector, need to subtract 1 from everything because bleaching classes start from 0
## sector = North:
##  habitat mean.class     SE  df asymp.LCL asymp.UCL
##  C             4.63 0.0467 Inf      4.54      4.72
##  F             4.52 0.0702 Inf      4.39      4.66
##  L             4.15 0.0846 Inf      3.99      4.32
##  U             4.90 0.0455 Inf      4.81      4.99
## 
## sector = Central:
##  habitat mean.class     SE  df asymp.LCL asymp.UCL
##  C             3.88 0.1186 Inf      3.65      4.11
##  F             3.48 0.1385 Inf      3.21      3.75
##  L             3.44 0.1470 Inf      3.16      3.73
##  U             4.06 0.1933 Inf      3.68      4.44
## 
## sector = South:
##  habitat mean.class     SE  df asymp.LCL asymp.UCL
##  C             2.11 0.0556 Inf      2.00      2.22
##  F             1.90 0.0615 Inf      1.78      2.02
##  L             1.56 0.0647 Inf      1.44      1.69
##  U             2.37 0.0938 Inf      2.19      2.56
## 
## Confidence level used: 0.95
emmeans(hughes.clmm, ~habitat|sector, mode = 'mean.class') %>%
  pairs() #to directly compare the habitats by sector
## sector = North:
##  contrast estimate     SE  df z.ratio p.value
##  C - F       0.105 0.0571 Inf   1.838  0.2557
##  C - L       0.476 0.0675 Inf   7.063  <.0001
##  C - U      -0.269 0.0555 Inf  -4.849  <.0001
##  F - L       0.372 0.0823 Inf   4.513  <.0001
##  F - U      -0.374 0.0748 Inf  -5.005  <.0001
##  L - U      -0.746 0.0881 Inf  -8.464  <.0001
## 
## sector = Central:
##  contrast estimate     SE  df z.ratio p.value
##  C - F       0.397 0.1062 Inf   3.740  0.0011
##  C - L       0.435 0.1076 Inf   4.045  0.0003
##  C - U      -0.180 0.1695 Inf  -1.061  0.7134
##  F - L       0.038 0.1360 Inf   0.279  0.9924
##  F - U      -0.577 0.1890 Inf  -3.053  0.0122
##  L - U      -0.615 0.1916 Inf  -3.210  0.0073
## 
## sector = South:
##  contrast estimate     SE  df z.ratio p.value
##  C - F       0.212 0.0457 Inf   4.641  <.0001
##  C - L       0.545 0.0488 Inf  11.173  <.0001
##  C - U      -0.265 0.0801 Inf  -3.313  0.0051
##  F - L       0.333 0.0557 Inf   5.981  <.0001
##  F - U      -0.477 0.0860 Inf  -5.547  <.0001
##  L - U      -0.810 0.0861 Inf  -9.412  <.0001
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
newdata <- emmeans(hughes.clmm, ~habitat|sector, 
                   mode = 'mean.class') %>%
  as.data.frame() %>%
  mutate(across(c(mean.class, asymp.LCL, asymp.UCL), function(x) x-1))

ggplot(newdata) +
  geom_hline(yintercept = 1, linetype = 'dashed', size = 0.1) +
  geom_hline(yintercept = 2, linetype = 'dashed', size = 0.1) +
  geom_hline(yintercept = 3, linetype = 'dashed', size = 0.1) +
  geom_pointrange(aes(y = mean.class, x = habitat,
                      ymin = asymp.LCL, ymax = asymp.UCL)) +
  facet_grid(~sector) +
  scale_y_continuous('Bleaching score', breaks = (0:4), labels = 0:4,
                     limits = c(0, 4), expand = c(0, 0)) +
  theme_bw()